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We present and test a new algorithm for time-evolving quantum many-body systems initially 
proposed by Holzner et al. [Phys. Rev. B 83, 195115 (2011)]. The approach is based on merging the 
matrix product state (MPS) formalism with the method of expanding the time-evolution operator 
in Chebyshev polynomials. We calculate time-dependent observables of a system of hardcore bosons 
quenched under the Bose-Hubbard Hamiltonian on a one-dimensional lattice. We compare the new 
algorithm to more standard methods using the MPS architecture. We find that the Chebyshev 
method gives numerically exact results for small times. However, the reachable times are smaller 
than the ones obtained with the other state-of-the-art methods. We further extend the new method 
using a spectral-decomposition-based projective scheme that utilizes an effective bandwidth signif¬ 
icantly smaller than the full bandwidth, leading to longer evolution times than the non-projective 
method and more efficient information storage, data compression, and less computational effort. 


I. INTRODUCTION 

Besides being of central interest in the field, achieving 
large accessible times in the time evolution of strongly- 
correlated quantum many-body systems with existing nu¬ 
merical methods has proven to be a daunting task in any 
spatial dimension and particularly for global quenches. 
Experiments in quantum many-body physics in the last 
years have evolved in such a manner that local control 
over degrees of freedom has become more feasible^^® and 
in which quantum magnetism, spin dynamics, and relax¬ 
ation dynamics have been explored. Additionally, along 
this experimental work a whole body of theoretical inves¬ 
tigations has arisen that relies on various analytical and 
numerical methods to describe the dynamics of these ex¬ 
periments. 

The time-dependent Schrbdinger equation 

( 1 ) 

for a generic time-independent Hamiltonian H and ini¬ 
tial state I'f/'o) = IV'(O)) is formally solved by the time 
evolution operator 

U[t) = exp(—iJTt), (2) 

where the reduced Planck constant h is set to 1. The 
time-evolved quantum state for arbitrary times is then 
given by 

IV'W) = = exp(-iiJt)|^/:o). (3) 

However, for a many-body quantum system, the dimen¬ 
sion of the Hilbert space grows exponentially with the 
number of constituents in the system under considera¬ 
tion, making it impossible to calculate the matrix ex¬ 
ponential in Eq. (3) exactly and, therefore, approximate 
methods are required. 


The purpose of this work is to discuss achievable evo¬ 
lution times for complex quantum many-body systems 
such as global quenches relevant to the current exper¬ 
imental efforts in the field. As the errors encountered 
in experiments are usually larger than those in numeri¬ 
cal calculations, we are not interested in an increase in 
accuracy. 

One method that has proven extremely useful is 
which is based on the description of 
the quantum state in terms of matrix product states 

(MPS)13-19 

(4) 

{'^1 {c^l 

where a = {cti ... ctat} is the computational basis, 
and are U-dimensional row and column vectors, re¬ 
spectively, and (i = 2,..., N — 1) is a D x D matrix. 
Theoretically, every quantum state can be represented by 
an MPS if infinite matrix dimensions are allowed^°. The 
practical relevance of such a state description lies in the 
fact that one can often very well approximate the exact 
quantum state by an MPS with finite matrix dimension. 
From this perspective, MPS presents a class of states that 
compress exact many-body quantum states such that the 
number of coefficients needed to describe the state scales 
linearly in the number of constituents as opposed to the 
exponential scaling in the exact representation. Further¬ 
more, the approximation made in the compression step 
is well understood^^ and can be controlled by the matrix 
dimension D. 

With the help of MPS, several methods have been 
developed to calculate the time evolution of one¬ 
dimensional many-body quantum systems^®. The ear¬ 
liest methods utilize the Trotter^^’^^ decomposition of 
the time-evolution operator. Later approaches approxi¬ 
mate the matrix exponential in the Krylov^^ subspace. 
Both methods have been successfully applied to a series 
of different physical problems. 
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Nevertheless, the times reachable with current meth¬ 
ods are still very limited making the development of new 
methods still a very important endeavor. The limitation 
of evolution times accessible with MPS-based methods 
is closely related to the amount of entanglement in the 
quantum state. The maximal entanglement between two 
subsystems describable by an MPS is given by the log¬ 
arithm of the matrix dimension D. On the other hand, 
it has been shown that the entanglement after a quan¬ 
tum quench grows typically linearly in time^^ leading to 
an exponentially-growing matrix dimension, which is re¬ 
quired in order to keep the error fixed. 

In this work, we test a new method for calculating the 
time evolution of one-dimensional quantum many-body 
systems as it was proposed in Ref. 25 by Holzner et al. 
We attempt to merge MPS with the method of approx¬ 
imating the time-evolution operator in terms of Cheby- 
shev polynomials. The procedure of expanding the time- 
evolution operator in terms of Chebyshev polynomials is 
general and requires in principle solely a matrix-vector 
multiplication. The MPS approach together with the 
representation of the Hamiltonian as a matrix product 
operator provides an efficient way to perform these oper¬ 
ations in the quantum many-body framework. A related 
approach based on Chebyshev polynomials has recently 
been successfully applied in the frequency domain to ob¬ 
tain an efficient impurity solver for the dynamical mean- 
field theory (DMFT) algorithm^®^^® and for calculating 
spectral functions^® and Green’s functions^®. 

For real-time dynamics in MPS, however, the approach 
has not been tested so far. In this paper, we test the 
new method (dubbed f-CheMPS) for a non-trivial sys¬ 
tem of hardcore bosons which evolve in time under the 
Bose-Hubbard Hamiltonian in a one-dimensional lattice. 
We show that time-dependent observables can be calcu¬ 
lated numerically exactly with the t-CheMPS method up 
to a certain time beyond which exponentially growing 
errors become dominant. The time reachable is given 
by the amount of entanglement in the n-th Chebyshev 
vector and can be slightly increased by making use of 
a projection procedure onto the energy range where the 
initial state has finite nonzero spectral weight. We com¬ 
pare our results to the time-evolution methods based on 
the Trotter^^’^^ expansion and Krylov^^ approximation 
of the time-evolution operator. We find that for the prob¬ 
lem considered in this work the Trotter-based method 
reaches the longest times, followed by the method based 
on the Krylov approximation. 

This paper is structured as follows: Section H gives a 
brief overview of the standard state-of-the-art methods in 
time evolution within the MPS context. Section HI dis¬ 
cusses the t-CheMPS method and its workings. Section 
IV presents an extension of the latter, namely, projective 
t-CheMPS based on the spectral decomposition of the ini¬ 
tial state. Section V discusses the Bose-Hubbard-model 
global quench used for the simulations in this paper, the 
results of which are documented in Section VI. The pa¬ 
per concludes with Section VH. 


II. STANDARD TIME-EVOLUTION METHODS 
IN MPS 

A. Krylov time-evolution 

Instead of treating Schrodinger’s equation as a differen¬ 
tial equation, one considers, for time-independent Hamil¬ 
tonians, the time-evolution operator exp(—iiJf). This 
sets the nontrivial task of evaluating an exponential of 
matrices^^’^^. One of the most efficient methods is the so- 
called Krylov subspace approximation^^d6,23^ where one 
realizes that our interest lies in eyip{—iHt)\'ip) rather than 
exTp{—iHt). In DMRG H\ip) is available efficiently, and 
this can be utilized through forming the Krylov subspace 
by successive Gram-Schmidt orthonormalization of the 
set {IV’), —■ ■ ■ }) where |^) is assumed 
to be normalized here. Here, —\Ht is approximated 
regarding its extreme eigenvalues by YTY"^, where Y 
is the matrix containing the n Krylov vectors thus ob¬ 
tained from the Gram-Schmidt decomposition and T is 
an n X n tridiagonal matrix. This approximation is up to 
a very good precision even for relatively small numbers of 
Krylov vectors^^. Thereafter, the exponential is given by 
the first column of Y expT, where the latter exponential 
is now much easier to calculate. 


B. Suzuki-Trotter time-evolution 

Another prominent and very efficient method for eval¬ 
uating the above matrix exponential is the (Suzuki- 
)Trotter decomposition^^’^®’^^’^^. This method is mainly 
useful for Hamiltonians with nearest-neighbor interac¬ 
tions. In the case of a one-dimensional chain, the Hamil¬ 
tonian H = Hi + H 2 is divided into odd- and even-bond 
terms, Hi and H 2 -, respectively, where Hi = J2^i ^ 2 i-i 
and H 2 = Here, hi is the local Hamiltonian 

linking sites i and i + 1^ and N is the total number of sites 
on the lattice. ^ 0 as neighboring local Hamil¬ 

tonians do not commute in general, but all the terms in 
Hi and H 2 commute. As such, the first-order Trotter de¬ 
composition of the infinitesimal time-evolution operator 
is 

Moreover, the second-order Trotter decomposition reads 

g-iilAt ^ g-iAiAt/2g-iA2Atg-iAiAt/2 (g) 

One can go for yet higher orders and conclude that an 
n*^-order Trotter decomposition will yield over a time 
step At an error of the order of (At)"+^. As one requires 
t/At time steps in order to reach an evolution time t, the 
error grows at worst linearly^^ in time t, and therefore, 
the resulting error is bound by an expression of the order 
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of (Af)"'f. For the purposes of this study, it turns out 
that second-order Trotter decomposition is optimal. 

Time-dependent DMRG (f-DMRG) uses adaptive 
Hilbert spaces that follow the state being optimally 
approximated, and was first proposed independently in 
the works of Daley, Kollath, Schollwock, and Vidal^^ and 
White and Feiguin®, based on the time-evolving block- 
decimation (TEBD) algorithm^°’^^ for the classical sim¬ 
ulation of the time evolution of weakly-entangled quan¬ 
tum states. Shortly afterwards, Schmitteckert^"^ pub¬ 
lished on nonequilibrium electron transport in interacting 
one-dimensional spinless Fermi systems using f-DMRG. 


III. t-CheMPS 



-I—i - ^ ^ 

-1 " 2 ^“^ " 1 


In this section, we review a recipe for time evolu¬ 
tion using the Chebyshev matrix product state approach, 
namely t-CheMPS, as it was proposed in Ref. 25 by 
Holzner et al. As such, a brief expose on Chebyshev 
polynomials is in order. 

Chebyshev polynomials of the hrst kind, Tn{x); n G N 
are given by the recursive relations 


FIG. 1. (Color online) (a) The spectral decomposition S{uj) 
of I'i/ifl) relative to H has nonzero weight in the region [E*,E*] 
where the effective bandwidth W* = E* — E* is significantly 
smaller than the full many-body bandwidth W — Es — Eg. 
(b) In the Chebyshev expansion approach for time evolution, 
it may be advantageous to rescale H by mapping the effective 
bandwidth from [E;,E*] to [-W',W'] where W' = 1 - f 
with et = 0.025 being a safety factor^®. 


f 1 

for n = 0, 

Tnix) = { X 

for n = 1, (7) 

[ 2xTn-l{x) - Tn- 2 ix) 

for n > 1. 


A useful non-recursive expression for the Chebyshev 
polynomials is 


T„(x) = cos(n arccosa;). (8) 

Moreover, they form an orthonormal set of polynomials 
on the interval x € [—1,1] with respect to the weighted 
scalar product 


dx 

(T„,r™)=/ ^ - M x)Trn{x), (9) 

J~i TTV1 — a; 

and are divergent in the region x G (—oo, —1) U (1, oo). 

Chebyshev polynomials have been extensively studied 
in the mathematics and engineering literature^®^^®. 


Since the Chebyshev polynomials of first kind are di¬ 
vergent outside of the region x G [—1,1], and knowing 
that these polynomials will be functions of a Hamilto¬ 
nian, this effective bandwidth is rescaled to [—IF', IF'] 
where IF' = 1 — ^ safety factor to guaran¬ 

tee that the domain of the Chebyshev polynomials will 
remain within I = [—1,1]. In our numerical simulations, 
tt has been set to 0.025. This rescaling, when applied to 
the original Hamiltonian H will lead to a rescaled Hamil¬ 
tonian H' where 


H' = 


H-b 

a 


( 10 ) 


with a = W* 1(2 — et) and b = [E* + E*)/2. Now one 
can express the Chebyshev polynomials of the first kind 
in terms of this rescaled Hamiltonian. 

Several constructions of Chebyshev approximations^® 
can be found for a given function /(a;)|xGi, with the one 
most suited for our purposes being 


A. Chebyshev expansion in the time domain 

We consider a system in which a Hamiltonian H acts 
on an initial state l'0o), thus propagating its time evo¬ 
lution. The full many-body bandwidth of IT is IF = 
Eg — Eg, where Eg [Eg) is the groundstate (skystate) en¬ 
ergy of IT. In many cases, this bandwidth is far larger 
than the effective bandwidth IF* = E* — E* that one can 
determine from the spectral function of |'0o) relative to 
H. As illustrated in Fig. 1, the spectral function of |'0o) 
relative to H has nonzero weight mainly over [E*,E*]. 


fix) 


1 

7r\/l — X® 


OO 

Mo ^ ^ , 

n—l 


where the Chebyshev moments Mn given by 


( 11 ) 


En = J f{x)Tnix)dx. (12) 

An order of N approximation / 7 v(a;) of f(x) is possible 
if one has access to the first N terms (0 < n < A — I), 
and thus it follows that 
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fN{x) 


1 

7r\/l — a;^ 


N-1 

n—l 


(13) 


The time-evolution operator can be expressed as {h = 1) 


This procedure of recurrence fitting effects variational 
minimization through a sequence of sweeps back and 
forth along the chain that proceed until the state be¬ 
ing optimized becomes stationary. Calling the state |t„) 
and after and before a fitting sweep, it becomes sta¬ 
tionary once the term 


U{t) = = J - II'), 

(14) 

and upon expressing the i5-function term therein as per 
Eq. (13), one obtains 

Af-l 

UN{t) = ^ (15) 

n=0 


Ac 


(ira l^ri) 

|lra) I I’ll l^ra) 


(19) 


drops below a specified fitting convergence threshold?^, 
which we have determined to suffice when set to 10“® for 
our purposes. 


C. Energy truncation 


with <f>o{t) = co(t) and <pn>o{t) = 2c„(t), where 

/ I -iatui'rp / i\ 

- J=^^dw' = (16) 

-1 TTV 1 — 

and Jn{at) is the Bessel function of the first kind of order 

n. 

It is to be noted here that in t-CheMPS, as will be 
elucidated later, one does not have to calculate the ac¬ 
tual wavefunction |^Ar(t)) = UN{t)\'fo) at a Chebyshev 
order N in order to determine the time evolution of some 
observable. 

B. Recipe for time evolution of initial state |i/)o) 

Here, we provide the steps needed to time-evolve an 
initial state |^/’o) under a Hamiltonian H using the t- 
CheMPS method. As an initialization step, we calculate 
the groundstate Ig) and the skystate |s) of H, noting 
that the skystate of H is nothing but the groundstate of 
—H. This allows us to determine the bandwidth W of H, 
from which we can make a specific choice for W* using 
the spectral-decomposition technique highlighted in the 
next section. Then we can determine a and b and rescale 
H to H' as per Eq. (10). The first Chebyshev vector |to) 
is set to the initial state |V'o)j while the second Chebyshev 
vector is given by \ti) = II'\to). Thereon, any Chebyshev 
vector |tn> 2 ) is obtained via the recursive relation 


The DMRG truncation step in the recursive-fitting 
procedure where H' is applied onto in order to cal¬ 

culate \tn) (see Eq. (17)) is not performed in the eigenba- 
sis of H', and, as such, high-energy components can be 
possibly passed on to subsequent recursion steps, lead¬ 
ing to divergences in higher-order Chebyshev vectors^^. 
This is remedied via energy truncation sweeps that oc¬ 
cur locally at each site through building the correspond¬ 
ing Krylov subspace, proceeding with the energy trun¬ 
cation at the site, and completing it before moving on 
to the next site. As DMRG truncation occurs in the 
recurrence-fitting procedure, no such further truncation 
is carried out here. The energy eigenbasis of H', where 
the energy truncation is to be performed, is not possible 
to access in full, and thus a Krylov subspace of dimension 
dx is constructed at each site. Then, a method such as 
Arnoldi’s algorithm is utilized to calculate the extreme 
eigenvalues of H' that are bigger than an energy trun¬ 
cation error bound per time step e in magnitude, where 
one can set e = 1.0, and focus on the proper value of W* 
based on the spectral decomposition of the initial state 
|'0o) with respect to II'. This is due to the fact that 
whatever value of W* one picks, the range of effective 
eigenenergies [A*, A*] will be rescaled to [—W, W], and 
in the Chebyshev context, the maximum and minimum 
energies must be no larger than e in magnitude. Further 
details on this method can be found in Ref. 25. 


D. Computing the time evolution of an observable 


|t„) = 2ff'|t„_i) - |t„_2). (17) 

This recurrence relation can be implemented using the 
compression or fitting procedure^®’^®. This procedure 
finds an MPS representation for |t„) by variationally min¬ 
imizing the fitting error"^^ 

Afit = |||t„) - (2ff'|t„_i) - |tn-2))||". (18) 


Consider that we wish to compute the time evolution 

= {il){t)\dj\il){f)) (20) 

of some observable O at a given site j on the chain. 
We represent the time-evolved state |^/’(t)) in terms of 
the Chebyshev representation of order N of the time- 
evolution operator of Eq. (15) on the initial state |'0o): 
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N-l 

^ r„(i^')</>.(i)IV'o), (21) 

n—0 

Noticing that IV'o) = l^o) and that Tn{H')\to) = |t„), Eq. 
(21) becomes 


decomposition of the initial state |'0o): and of its time- 
evolved version, Itpit)). We now discuss how this can be 
done, focussing first on |'0o)j and thereafter generalizing 
the discussion to |V'(i)) in the Appendix. For ease of no¬ 
tation, in this Section and the Appendix, H shall denote 
the rescaled Hamiltonian of our system with bandwidth 
W = Eg — Eg, and w £ [—1,1]. 


N-l 

^ ijinmn). (22) 

n—0 

Plugging Eq. (22) into Eq. (20), we get 

N-l 

{Oj){t)= (23) 

n,m—0 

As already mentioned, in t-CheMPS one is never ob¬ 
ligated to calculate the actual wavefunction lipit)) itself 
in order to calculate a certain observable using Eq. (23). 
Furthermore, the coefhcents cj)n{t) = {—i)^Jn{at) {n > 0) 
decay rapidly with n for n > at. It is therefore possible 
to define a maximaum time for a given number of Cheby- 
shev moments N such that the neglected weight in terms 
of the coefficients is smaller than a certain threshold. We 
define tmax as the largest t such that 


^ K{t)<t>n{t) < 10-3, (24) 

n=N 

which is justified because the moments (tm\Oj\tn) decay 
quickly with |n — m| [see Fig. 8]. In practice we deter¬ 
mine tmax by calculating J2n=N 4>n(t)4’n{t) < I0“3, with 
A^max = 500 for which we have (/'Ar„„(t) < 10-3°° in 
the relevant time range or (f’Nma.N) = 0 all practical 
purposes. 


IV. PROJECTIVE t-CheMPS 

We wish to find a way to calculate the effective band¬ 
width of a wavefunction ['0 (t)) at a time t. The reason 
behind this is that in t-DMRG one uses the full band¬ 
width while time-evolving the wavefunction and that 
leads to smaller evolution times that can be reached nu¬ 
merically. Using a smaller effective bandwidth may lead 
to larger numerically-accessible evolution times. 

Let the full many-body bandwidth of the model be 
W = Eg — Eg. Suppose that the initial state |0o) = 
|0(t = 0)) has spectral support on a limited frequency in¬ 
terval [Eg, Eg], of width W* = E* — E*, where E* < Eg 
and E* > Eg. Then, it would be possible to do the time 
evolution with t-CheMPS by rescaling this effective band¬ 
width, rather than the full bandwidth, onto the interval 
[—1,1]. Thus, it is of interest to explore the spectral 


Spectral decomposition of initial state |0o) 
The spectral decomposition of |0o) is 


'S'(w) = (0o|<5(w - iL)|0o)- (25) 

A Chebychev expansion of the (5-function of order N has 
the form: 


5n{i^ - H) = 


r\/l — I 


7V-1 


go + 2 ^ gnTn{H)Tn{uj) 


, (26) 


where the coefficient g„ is a Jackson damping coefficient 
defined as 


9n = 


{N — n + 1) cos -I- sin 


cot ^ 
Af-l-l N+l 


N+\ 

We introduce 0„ such that 


(27) 


f go if n = 0 , 

\ 2g„ if n > 0. 
This allows us to write Eq. (26) as 


(28) 


JV-l 

Sn{^ -H)= . -^ V 0nTn{H)T^{io). (29) 

1 - n=0 

Now we calculate S{uj) using Eq. (29) and noting that 
our initial wavefunction |'0o) equals the first Chebyshev 
vector [to) and that |t„) = T„(iL)|to): 


S{uj) = (0o|'5Ar(w - Hfifo) 

N-l 

= -2 E ^«ir„((u)(to|t„), (30) 

1 - 


Hence, all we have to do to calculate S'(w) is to cal¬ 
culate the moments \i.n = {to\tn). To achieve a specified 
spectral resolution of, say. A, we have to use an expansion 
order of Va = 2IT/A. Moreover, we provide in the Ap¬ 
pendix a derivation in terms of the Chebyshev moments 
of the spectral decomposition of the time-evolved wave- 
function \'tp{t > 0)), which can be used as a numerical- 
fidelity check. 
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2 3 4 5 

time t (1/J) 


FIG. 2. (Color online) The time evolution for the particle 
density at site Z//2 after a global quench with (a) [7 = 0, 
(b) [7 = 2, and (c) U = 5, each obtained with the Trot¬ 
ter (blue dotted), Krylov (green dashed) and t-CheMPS (red 
solid) methods. For [7 = 0 the exact time evolution given 
by Eq. (36) is shown in light grey. All the methods give nu¬ 
merically exact results for short times. The time reached by 
the t-CheMPS method is given by tmax (see Eq. (24)). After 
this time threshold the error quickly increases. In all cases 
the Trotter method reaches the longest times followed by the 
Krylov method. 


L/2 

|^o) = |^(0)) = n^Li|0) (31) 

i=l 

that is a bosonic lattice of size A = 32 in which every 
odd site has a single boson and every even site holds 
zero occupancy. It can be thought of as the groundstate 
of some suitable Hamiltonian. The system is globally 
quenched to the Bose-Hubbard Hamiltonian 

L — 1 jj L 

H = -J^ {h%+i + h.c)j + ^ X! 

i=l ^ i=l 

where J and U are the hopping and interaction terms of 
the Bose-Hubbard model. This global quench has already 
been studied using . We consider different 

values of the on-site interaction strength, including the 
analytically solvable case of [7 = 0. 

At [7 = 0, the Hamiltonian in Eq. (32) reduces to 

L-l 

H = - + h.c.^ . (33) 

In the case of non-interacting bosons {U = 0), scat¬ 
tering is not the physical mechanism behind local relax¬ 
ation. Instead, the time-dependent contributions to the 
reduced density operator of the regarded subsystem con¬ 
sist of quickly oscillating phases that average out under 
sufhcient conditions leading to a relaxation of the density 
operator^^. In the current case, excitations start propa¬ 
gating from all sites with a finite speed throughout the 
duration of the time evolution spreading the information 
about the initial conditions more and more over the en¬ 
tire system. The incommensurate mixing of these excita¬ 
tions then can lead to a state that appears to be locally 
perfectly relaxed. This case leads to an exact analytical 
solution covered in Ref. 42 by a Fourier transformation of 
the ladder operators involved. In the Heisenberg picture 
the time evolution of the ladder operators reads 


V. GLOBAL-QUENCH TEST MODEL 


For the comparison we wish to carry out between the 
Suzuki-Trotter decomposition, the Krylov approximation 
and the t-CheMPS methods, we consider a benchmark 
test model: a strong global quench in the Bose-Hubbard 
model (BHM) on a bosonic lattice at half filling with 
odd-site unity filling for different values of the on-site 
interaction strength U. Global quenches happen when 
an initial state undergoes a time evolution due to a new 
Hamiltonian for which the initial state has an extensively 
different energy as for the original Hamiltonian whose 
groundstate it was. 

We consider an initial state 


h{t) = 6,(0), 


k 1^1 
L 




k{l—i) —2\J cos 


wSKo), 


k 1^1 


(34) 

(35) 


where k = ^l, where I = 1,2,...,L. As ni{t) = 
one obtains 


l(l + (-l)*+ij„(4Jt)). (36) 
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FIG. 3. (Color online) Density-density correlations ^j{t) = for [7 = 2 (black) and U = 5 (blue) 

are shown for different half distances j: (a) j = 3, (b) j = 5, (c) j = 7, and (d) j = 9. All the correlators are obtained 
by the Trotter (dotted), Krylov (dashed), and t-CheMPS (solid) methods. All the methods give the same results up to the 
corresponding reachable times. 


5 ? 



FIG. 4. (Color online) The spectral decomposition S{u!) of the 
initial state I'tpo) as defined in Eq. (30) obtained by the proce¬ 
dure explained in Sec. IV. The green solid line shows the spec¬ 
tral decomposition of |^o) at interaction strength U = 2, the 
dashed blue line shows the spectral decomposition at [/ = 5. 
In both cases the spectral weight is located at the lower end 
of the spectrum and is negligible at higher energies. This al¬ 
lows one to determine proper projections for reduced effective 
bandwidths and to use the projective t-CheMPS method. 


VI. RESULTS AND DISCUSSION 

A. Convergence and performance 

A first quantity to gauge the convergence parameters of 
all three methods under consideration is provided by the 
particle density of this global quench with U = 0. The 


results shown in Fig. 2 exhibit good convergence for our 
purposes where the maximum on-site occupation num¬ 
ber in t-DMRG is set to (n)niax = 10 in accordance with 
Ref. 42. It is worth mentioning at this point that the 
underlying intention behind this work is not to contrive 
a method that surpasses standard methods such as Trot¬ 
ter time evolution and Krylov time evolution in terms 
of accuracy, as the latter have proven to be very precise 
with the right set of parameters in place. However, the 
goal is to investigate whether an alternative method such 
as t-CheMPS can, at the same accuracy or that within 
what is acceptable from an experimentally-suitable point 
of view, achieve larger times than those possible in Trot¬ 
ter time evolution or Krylov time evolution, especially 
that it has been demonstrated that Chebyshev polyno¬ 
mials can be very useful in time evolution at least outside 
of the context of MPS^^. 

In this work, we use a 2”‘^-order Trotter decomposition 
as in previous work^^ it has shown to be far more efficient 
than either 1®*- or 4*^-order Trotter decompositions in 
terms of accuracy and computational effort, respectively, 
while achieving approximately the same evolution times. 
The Krylov method we use employs an Arnold! iteration, 
which is considered to be the most efficient in Krylov 
implementations^®. In our calculations, we find that 
our Trotter calculations are convergent for a time step 
At = 0.01/J and truncation or fidelity threshold^®’®^’"'® 
of 10“® for each time step, while Krylov and t-CheMPS 
calculations are convergent for a fidelity threshold of 10“® 
for each time step. As shown in Figs. 2 and 3 for the 
particle density and density-density correlations, respec¬ 
tively, Trotter decomposition is the best method when 
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time t (1/J) time t (1/J) 


FIG. 5. (Color online) Projective t-CheMPS results for the particle density at L/2 of the global quench with (a) t/ = 2 and (c) 
U = 5 and different projection factors a = 0.5, 0.6, 0.7 for an effective bandwidth W* = aW as defined in Eq. (37), laid over 
the corresponding non-projective (a = 1) t-CheMPS results. Panels (b) and (d) show the same data as (a) and (c) but zoomed 
in to the relevant times, where errors start to diverge. Larger evolution times can be reached with the projected effective 
bandwidth in addition to less computational effort in the calculations and a smaller required disk space for representation of 
dynamics. 


it comes to largest accessible times. The times achieved 
by the Krylov method (shown in Fig. 2) match those 
arrived at by Flesch et al. in Ref. 42 for the same sys¬ 
tem. The accuracy of the t-CheMPS method is quite im¬ 
pressive and its results are actually quite exact for short 
times, but it can exceed neither the Krylov method nor 
the Trotter method in terms of largest evolution times 
reached. 


B. Projective t-CheMPS results 

As a further attempt at improving the results attained 
by the f-CheMPS method, we undertake the spectral de¬ 
composition in the cases of 17 = 2 ((n)max = 8) and 
U = 5 ((h)max = 4), the spectral functions of which are 
shown in Fig. 4. The case of 17 = 0 is not included as the 
bandwidth cannot be further reduced from its full size. 
Immediately, one notices that they both have nonzero 
weight mostly on the left half of the energy axis, i.e. in 
the lower energy half of the bandwidth. In the calcula¬ 
tions performed for this study, it has proven necessary 
to set E* = Eg to keep the Chebyshev approximation 
convergent, while E* G [Eg + W/2,Es). Here, Es is pro¬ 
jected according to 

E; = Eg + a- W, (37) 

where a € (0,1] is the projection factor, and when a = 1, 
it is in fact non-projective l-CheMPS that is being used 


and energy truncation is turned off in the simulations. 

Reducing the full bandwidth of the system to a re¬ 
duced effective bandwidth may lead to less computa¬ 
tional effort as one now requires fewer Chebyshev vec¬ 
tors in order to reach a certain maximum evolution time 
fmax, which is related to the expansion order ATnax by 
approximately^® tmax ~ A^nax/o, though, for our pur¬ 
poses, it is slightly smaller (see Eq. (24)) in order to 
achieve a desired precision (Sec. HID). Since a scales 
proportionally with the reduction in the full bandwidth 
upon projection, one need only achieve the same number 
of vectors in projective f-CheMPS as in non-projective t- 
CheMPS to facilitate a maximum evolution time bigger 
by that same factor of reduction in the full bandwidth. 
However, in projective f-CheMPS a new function enters 
into the computation, namely that of energy truncation. 
In our numerical simulations, the projective f-CheMPS 
method at any factor of reduction is unable to calculate 
up to the same order of expansion as the non-projective 
t-CheMPS method due to the additional computational 
effort of energy truncation, but, nevertheless, for certain 
projections, evolution times bigger than those achieved 
in non-projective f-CheMPS are reached as shown in Fig. 
5 for both cases U = 2 and U = 5. The best result is 
attained for both U values at a = 0.5, the most strin¬ 
gent projection factor used that did not lead to diver¬ 
gences, where an improvement of 20% (12%) is achieved 
for ?7 = 2 (C/ = 5) in terms of largest accessible evolution 
times. The corresponding converged Trotter results are 
overlaid for reference. It can be seen that even though 
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projective t-CheMPS does indeed reach greater evolution 
times than its non-projective counterpart, it still does 
not improve over the Trotter or Krylov methods. It is 
also worth noting here that projective t-CheMPS, despite 
even sometimes significant reductions in the full band¬ 
width of the system, still offers exact results for short 
times. 

Though one may be tempted to think that the pro¬ 
jective t-CheMPS method must achieve longer times the 
more one projects (be., the smaller a is), this is not the 
case in reality as then more computational effort is re¬ 
quired by the energy-truncation module that at some 
point it simply cannot handle all the required energy pro¬ 
jections when a <C 1, and in fact this renders the max¬ 
imum evolution time reachable smaller than that in the 
non-projective t-CheMPS method or it may outright lead 
to divergences^®. On the other extreme, if a < 1, then 
W* < W, and thus the computational effort is almost the 
same as in the non-projective t-CheMPS method with the 
added cost of energy truncation, which leads to evolu¬ 
tion times shorter than those attained by non-projective 
t-CheMPS. Thus, one has to choose a in a manner where 
energy truncation is not pushed to its limits and while at 
the same time W* is nontrivially smaller than W. 


C. Middle-bond dimension and vector size 

To avoid any confusion, we remind the reader here that 
the expansion order N is the total number of Chebyshev 
vectors, and each of the latter is indicated by an expan¬ 
sion index n that goes from 0 for the first vector to iV — 1 
for the highest-index vector. 

Intuitively, the projective t-CheMPS vectors are ex¬ 
pected to comprise of higher bond dimensions at the same 
expansion order than their non-projective t-CheMPS 
counterparts as exhibited in Fig. 6(a) and (b) for C7 = 2 
and U = 5, respectively, at the middle or central bond. 
This can be attributed to the fact that for the same time 
t attained in both methods, expansion order N* achieved 
by projective t-CheMPS for faithful representation of the 
system dynamics at this time is related to the corre¬ 
sponding expansion order N attained by non-projective 
t-CheMPS through N* m a ■ N < N. Hence, the vec¬ 
tor of a certain expansion index carries more information 
when generated by projective rather than non-projective 
t-CheMPS, because the generated N* Chebyshev vec¬ 
tors in projective t-CheMPS still, assuming convergence, 
must carry the same information about the system as the 
N (Ri N*/a > N*) Chebyshev vectors in non-projective 
t-CheMPS do. However, it can be seen in Fig. 6(a) and 
(b) that the projective t-CheMPS vectors do not reach 
the maximum central-bond dimensions that occur for the 
non-projective t-CheMPS vectors. In principle, one may 
expect that all the t-CheMPS vectors, regardless of the 
value of a ought to reach the same maximal matrix di¬ 
mensions. This is only true, however, if the workings 
of these calculations are the same, but this is not the 


case because in projective {a < 1) t-CheMPS an addi¬ 
tional computation effort is needed, that of energy trun¬ 
cation, which does not occur in non-pro iective (a = 1) 
t-CheMPS. 

In addition to obtaining fewer vectors required to ar¬ 
rive at an evolution time t in projective t-CheMPS, one 
finds that these vectors are in fact smaller in size the big¬ 
ger the reduction in bandwidth, be., the smaller a, is. In 
Fig. 6(a) and (b), each temporal isoline is constructed for 
a properly selected evolution time t that is appropriately 
matched to its corresponding expansion orders Na for 
the different a values based on the criterion in Eq. (24) 
(one may equally well use the more relaxed criterion of 
Na ~ aWt/{2 — et), which our calculations show is also 
adequate for [/ = 2 and U = 5). These temporal iso¬ 
lines indicate that at a time t, the corresponding projec¬ 
tive and non-projective f-CheMPS maximum-expansion- 
index vectors |tAr*_i) and |tjv-i), respectively, are such 
that the latter has larger matrix dimensions than the for¬ 
mer, and this becomes more pronounced the larger t is. 
It is interesting to also look at the total sum of matrix 
dimensions at the central bond of the Chebyshev vectors 
involved in arriving at a time t in t-CheMPS for different 
values of a. If Dn = D(\tn)) indicates the matrix dimen¬ 
sion at the central bond of |t„), then ^ would 

be a good measure of the computational effort required 
to reach a time t based on Eq. (24) in (non-)projective 
t-CheMPS for some value of a. This measure not only 
incorporates the maximal matrix dimension attained by 
the highest-index Chebyshev vector required to reach an 
evolution time t, but it also accounts for how many vec¬ 
tors are required to reach t, and this number varies de¬ 
pending on the value of a. This measure is depicted in 
Fig. 6(c) and (d) for ?7 = 2 and U = 5, respectively, 
where one can conclude that the more reduced the effec¬ 
tive bandwidth is (the smaller a is), the smaller is the 
computational effort required to reach a certain evolu¬ 
tion time t. For the longest common time arrived at by 
all calculations (t ~ 2.3/J), there is a factor of roughly 4 
with regards to mitigation of computational effort from 
a = 1 to a = 0.5. 

Moreover, at a certain evolution time t corresponding 
to a set of expansion orders Na for the different a-valued 
t-CheMPS calculations, one finds that the highest-index 
Chebyshev vector \tN^-i) occupies less disk space the 
smaller a is. If dn = d{\tn)) is the disk space occupied by 
Chebyshev vector |f„), then J2n=Q^ total disk 

space needed to house those Chebyshev vectors required 
to arrive at the dynamics up to time t corresponding to 
Na as per Eq. (24). This is presented in Eig. 6(e) and 
(f) for U = 2 and U = 5, respectively, where it can be 
seen that the greater the projection (or the smaller a is), 
the more reduction one obtains in total disk space. In 
fact, at t = 2.3/J, the reduction is more than an order 
of magnitude from a = 1 to a = 0.5. Therefore, upon 
projection, one obtains fewer Chebyshev vectors that as 
a whole are also smaller in size while representing the 
same dynamics, which indicates data compression. 
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FIG. 6. (Color online) The matrix dimension at the central bond for the Chebyshev vectors due to the projective (a < 1) and 
non-projective (a = 1) t-CheMPS methods for U — 2 (left column) and [7 = 5 (right column). In (a) and (b), one notices that 
at any expansion index n, the corresponding Chebyshev vector \tn) carries a larger central-bond dimension the smaller a is 
(he., the more reduced the effective bandwidth is), as at the corresponding expansion order n-P 1, t-CheMPS represents longer 
dynamics the smaller a is. However, in (c) and (d), one notices that the sum of the matrix dimensions at the central bond 
for the Chebyshev vectors leading up to a certain time t is far smaller the lower the value of a, indicating less computational 
effort upon greater reduction in the bandwidth. Note how the projective t-CheMPS vectors cannot reach the maximum matrix 
dimension non-projective f-CheMPS vectors have, and this is due to the additional computational effort of energy truncation 
necessary in the projective t-CheMPS method but nonexistent in its non-projective counterpart. Additionally, (e) and (f) show 
significant conservation of disk space in projective t-CheMPS for any evolution time t, where the smaller a is, the less disk 
space is required for storing the Chebyshev vectors that are necessary to represent dynamics up to t. Note that the selected 
times are indicated via temporal isolines in (a) and (b). 


D. Comparison with other methods 

It is interesting to produce a quantitative comparison 
of the I-CheMPS method with the Krylov and Trotter 
methods in the time domain. One can consider the ma¬ 
trix dimension at the central bond required to faithfully 
represent the dynamics over the evolution times. One can 


again here represent the central-bond total matrix dimen¬ 
sion for t-CheMPS at a time t as J2n=o^ where is 
the matrix dimension at the central bond of |t„), and 
and t are related as per Eq. (24), as is done in Fig. 6, but 
this representation would not be fair as it encompasses 
all the Chebyshev vectors required to construct the wave- 
function \tjj[t)), the construction of which, unlike in the 
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FIG. 7. (Color online) The matrix dimension at the cen¬ 
tral bond of the bosonic chain using the t-CheMPS, Krylov, 
and Trotter methods for interaction strengths (a) t/ = 0, (b) 
U = 2, and (c) U = 5. The matrix dimension for t-CheMPS 
increases very quickly and greatly exceeds its counterparts in 
the Krylov and Trotter methods for common evolution times 
irrespective of the interaction strength. This prohibits the 
method from achieving much longer times. 


at those interaction strengths, while the non-projective 
t-CheMPS result is used for = 0 as there no projection 
is possible. The comparison is displayed in Fig. 7, where 
it can be noted that the central-bond matrix dimension 
required in the t-CheMPS method to represent the dy¬ 
namics up to a common evolution time t is much greater 
than that in the Krylov or Trotter methods regardless of 
what the interaction strength U is. This is in agreement 
with Ref. 29, particularly in the case where the rescaled 
Hamiltonian H' is simply a factor of the original Hamil¬ 
tonian H, which is the case in this study when U = 0, 
depicted in Fig. 7(a). At this interaction strength, the 
skystate and groundstate energies are equal in magni¬ 
tude but of opposite sign, rendering 6 = 0. This leads to 
H' = 77/a, which is the condition proven in Ref. 29 to 
assert that then the Chebyshev vectors are equivalent to 
time-evolved wavefunctions for a proper time step in the 
Krylov or Trotter methods^®. 


E. Matrix moments 

Finally, in Fig. 8, we take a closer look at the behavior 
of the Chebyshev moments {tm\nL/ 2 \tn)- In particular, 
we observe that these moments carry the greatest weight 
along the back diagonal (\). This is to be expected as 
these Chebyshev moments involve two states |t„) and 
\tm) that have little overlap, since, if m > n, \tm) is ar¬ 
rived at by consecutively applying H'm — n times onto 
|7„), and as the latter is not an eigenstate of 77', this ren¬ 
ders the two vectors with little overlap the bigger |m —n| 
is. Moreover, the cross sections of these moments along 
the dotted black lines in Fig. 8(a)-(d), corresponding to 
some expansion order, say N[a) = lOOa, exhibit a de¬ 
caying behavior around |m — n| =0 that is zero for large 
\m — n\. These cross sections for the different a values 
are depicted in Fig. 8(e). As mentioned previously, this 
validates the constraint for the maximum evolution time 
tmax that '/n(*)'/n(7) < 10“^ while neglecting off- 

diagonal terms as indeed one can see that the Chebyshev 
moments in Fig. 8 carry nontrivial weight mostly for 
quite small values of |m — n|, thereby making this con¬ 
straint sufficient. 


Trotter or Krylov methods, is never undertaken in the 
t-CheMPS method (see Sec. HI D) . Thus, even though 
this representation is proper in Fig. 6(c) and (d) as it 
involves a comparison between the different bandwidth 
reductions in the t-CheMPS method through covering 
the number of Chebyshev vectors involved in reaching 
an evolution time 7, for comparison with the Trotter and 
Krylov results, is the proper quantity to look at, 

because HjVq-i is the largest matrix dimension of the 
central bond attained by any of the Chebyshev vectors 
required for faithful representation of the dynamics up to 
evolution time t. For this comparison, the a = 0.5 projec¬ 
tive 7-CheMPS result for U = 2 and 7/ = 5 is chosen as it 
performs best compared to other 7-CheMPS approaches 


VII. CONCLUSION 

A new method, 7-CheMPS, based on the Chebyshev 
expansion in the time domain in the context of MPS 
has been presented for calculating the time evolution of 
quantum many-body systems, including global quenches. 
Using a test system of importance in the field of quan¬ 
tum many-body physics, we demonstrate that 7-CheMPS 
arrives at exact solutions of a given observable for short 
times, but does not exceed the largest times accessible 
by standard time-evolution methods such as the Trot¬ 
ter decomposition and the Krylov approximation. Fur¬ 
thermore, a projective version of the method, projective 
























12 



0.32 

0.24 

O.IG 

0.08 

0.00 

-0.08 

-0.16 

-0.24 


-0.32 


expansion order n 


expansion order n 



FIG. 8. (Color online) The behavior of the Chebyshev mo¬ 
ments {tm\nL/ 2 \t’n) for (a) non-projective (a = 1) and pro¬ 
jective t-CheMPS for (b) a = 0.7, (c) a = 0.6, and (d) 
a = 0.5. It can be seen how the bulk of the information 
lies where |m — n| < 15 and biggest around m « n. In (e) 
the cross sections of these moments along the back diagonal 
{{Na,0), {0, Na)} where N{a) = lOOo? are shown, displaying 
rapid decay around |m — n| = 0 for all a values. 


f-CheMPS, based on spectral decomposition and system- 
bandwidth reduction is introduced that improves on the 
largest evolution times accessible while significantly eas¬ 
ing computational effort and greatly reducing disk space 
for the same dynamics. Moreover, we find again that 
Trotter expansion is still the favorable method with re¬ 
gards to largest accessible evolution times. 
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Appendix: Spectral decomposition of the 
time-evolved wavefunction \ip{t > 0)) 

Reminding the reader that here for ease of notation, as 
in Sec. IV, H is taken to be the rescaled Hamiltonian with 
bandwidth W = Eg — Eg, and uj € [—1,1], we proceed 
with first remarking that the spectral decomposition of 
the time-evolved state |V'(i)) is equal to that of the initial 
state |'0o): 


{ip{t)\S{uj - H)\'ilj{t)) 

= (^o|yt(t)J(w-7L)U(t)|^o) 

= {iIjo\S{uj - H)\'ijjo), (A.l) 

where the time-evolution operator U(t) commutes with 
S{uj — H). Thus, it is a good check of the validity and 
convergence of the Chebyshev vectors to ascertain that 
the spectral decomposition is the same at any time t 
when calculated by the corresponding Chebyshev mo¬ 
ments. Furthermore, this may also be employed as an 
alternate way to Eq. (24) to determine how many Cheby¬ 
shev vectors one would need to faithfully represent the 
physics at an evolution time t, using the error with re¬ 
spect to S{uj) in Eq. (30) as a gauge. As such, we pro¬ 
vide here a derivation that allows one to calculate the 
spectral decomposition of the time-evolved wavefunction 
|'0(<)) from the corresponding Chebyshev moments. 

In the t-CheMPS method, one can represent the wave- 
function \tp{t)) = {/(Ql'i/'o) as 


N-l 

^ <^„(t)|t„) (A.2) 

n—0 

as per Eq. (21). To reach a specified evolution time t, we 
need to use an expansion order of Nt ~ tW 12 or, more 
stringently, as specified by Eq. (24). 

Now, to calculate the spectral function StiuS) for the 
wavefunction with a specified spectral resolution 

A, we can proceed as follows: 


Stiuj) = (?/>(t)|(57VA(w - H)\ip{t)) 


Nt-l 


= - H)\tn) 

n,n'—0 

Nt-l Na-1 

(A.3) 


where = (t„/|r„//(iJ)|t„). 

Note that we need to use different upper limits on the 
sums on n and n' than on the sum on n", in order to 
reach a specified time t with a specified spectral resolu¬ 
tion A. To evaluate the moments arising here, we recall 
the following identity: 
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(A.4) 

It is advisable to use it in such a way that the order of 
the polynomials that arise remain as small as possible. 
Thus, for the case that n < n', we proceed as follows 
(with ni = n" and n 2 = n): 

T„„{H)T„iH) = (A.5) 

which leads to 

= {tn^\Tn.{H)\tr,) = {tn'\Tn'>{H)Tn{H)\to) 

= + ■^{tn'\T\n-n"\iH)\to) 


For the case that n' < n, we proceed analogously, but 
with rii = n' and n 2 = n": 


T^,[H)T^.{H) = + \T\^'-n"\{H), (A.7) 

which in turn leads to 

= {tn’\Tn"{H)\tn) = {to\Tn'{H)T^n{H)\t^) 

= 2{'^o\Tn'+n"{H)\tn) + 2^to\T\n'-n" \{H)\tn) 

~ +n"\tn) + —n" \\in) ■ 

Thus, we conclude that in order to calculate St{ui) with 
a specified resolution of A up to a specified time t, we 
need all Chebyshev vectors up to order Aa + A^t- 
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